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Abstract 

We present a single step, second-order accurate Godunov scheme for ideal MHD 
which is an extension of the method described in [1] to three dimensions. This algo- 
rithm combines the corner transport upwind (CTU) method of Colella for multidi- 
mensional integration, and the constrained transport (CT) algorithm for preserving 
the divergence-free constraint on the magnetic field. We describe the calculation of 
the PPM interface states for 3D ideal MHD which must include multidimensional 
"MHD source terms" and naturally respect the balance implicit in these terms by 
the V • B = condition. We compare two different forms for the CTU integration 
algorithm which require either 6- or 12-solutions of the Riemann problem per cell 
per time-step, and present a detailed description of the 6-solve algorithm. Finally, 
we present solutions for test problems to demonstrate the accuracy and robustness 
of the algorithm. 
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1 Introduction 



In a previous paper [1] we described a two-dimensional, second-order accu- 
rate Godunov method for ideal MHD that evolves the magnetic field using 
the Constrained Transport (CT) [11] algorithm for preserving the divergence- 
free constraint on the magnetic field. In its simplest form, CT requires area- 
averaged values of the magnetic field which are stored at cell faces. We argued 
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that this is the most natural discrete representation of the field in that the 
integral form of the induction equation is based on area (rather than volume) 
averages, and therefore the discrete form of the equations should respect this 
difference. There are three important ingredients to our MHD algorithm: (1) 
a modification of the piecewise parabolic method (PPM) reconstruction step 
used to construct time-advanced estimates of the conserved variables on cell 
faces that are fed to the Riemann solver to incorporate multidimensional terms 
essential in MHD, (2) a new method for constructing the fluxes (at cell edges) 
of the area-averaged magnetic fields (at cell faces) from the fluxes returned by 
the Riemann solver (at cell faces) of volume averaged magnetic fields (at cell 
centers) which are based on the fundamental relationship between the area- 
and volume-averaged variables, and (3) a directionally unsplit integration al- 
gorithm based on the Corner Transport Upwind (CTU) method [7]. 

Through a series of test problems, we demonstrated the importance of each of 
the ingredients to our algorithm. In particular, we showed through tests based 
on the advection of two-dimensional field loops that our new methods for 
constructing the fluxes needed by the CT algorithm are essential for stability, 
and are an improvement over previous Godunov methods that use CT, e.g. 
[3]. We also showed that by using the second-order accurate CTU integration 
algorithm, a method could be constructed which has less numerical dissipation 
and has the important property of reducing exactly to the one-dimensional 
algorithm for plane-parallel, grid-aligned flows. Since CT does not require 
costly solutions to elliptic equations, we expect MHD Godunov schemes based 
on CT to be more cost effective that those that use divergence-cleaning [9,20]. 
Given the attractive properties of the method, it is of interest to extend it to 
three-dimensions for use in applications. 

When directional splitting is used, the extension of Godunov methods from 
two- to three-dimensions is usually trivial. However, directional splitting is 
unsuitable for MHD, because it is impossible to enforce the divergence-free 
constraint between partial updates unless all three components of the magnetic 
field are updated together, which in turn violates the assumption basic to 
splitting that each dimensional operator is independent and can be split from 
the others. As a result, in [1] we adopted the unsplit CTU integration scheme. 
Even in hydrodynamics, the extension of CTU to three-dimensions is not 
trivial [18]. For our MHD algorithm, extension to three-dimensions requires 
modifying two of the three ingredients of the method, in particular (1) the 
PPM reconstruction algorithm must be modified to include multidimensional 
terms for MHD in such a way as to respect a balance law implied by the 
V ■ B = condition, and (2) the CTU algorithm must be modified to include 
source terms as well as the transverse flux gradient terms. The primary purpose 
of this paper is to describe in detail these modifications and to demonstrate 
that the resulting algorithm is both accurate and robust. 
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We extend our MHD test suite to three- dimensions to demonstrate the accu- 
racy and fidelity of our method. We find that, once again, the passive advection 
of a multidimensional field loop is a challenging test of finite volume meth- 
ods for MHD. In particular, for a field loop confined to the (x, y)-plane in 
three-dimensions advected with a constant velocity with v z ^ 0, the vertical 
component of the magnetic field B z will evolve unless care is made to ensure 
the multidimensional balance of MHD source terms in both the PPM charac- 
teristic tracing step and the transverse flux gradient update step. In fact, this 
observation leads to a useful definition of the appropriate difference stencil on 
which the divergence-free constraint must be maintained. If V • B = on a 
stencil which is different from that used to construct the fluxes of B z , the latter 
will show unphysical evolution in this test for conservative algorithms. On the 
other hand, if a numerical method keeps B z constant to round-off error on the 
test, it must preserve the divergence- free constraint on the appropriate stencil. 
Moreover, this test is another demonstration that it is essential to maintain 
the divergence-free constraint exactly in MHD, as was originally emphasized 
by [6]. This test, along with linear wave convergence tests, a test based on the 
propagation and convergence of nonlinear, circularly polarized Alfven waves, 
and multidimensional blast wave tests are all presented in section 6. 

The paper is organized as follows. In section 2, we write down the equations of 
ideal MHD solved by our method, and describe the finite-volume discretiza- 
tion of mass, momentum, and energy, and the finite-area discretization of the 
magnetic field. In section 3, we describe our extension of the PPM reconstruc- 
tion algorithm to three-dimensional MHD. In section 4, we briefly review the 
upwind CT algorithms introduced in [1], while in section 5, we describe two 
formulations for the CTU integration algorithm to 3D. In section 6 we present 
the results of our test suite, while in section 7 we conclude. 



2 Ideal Magnetohydrodynamics and Constrained Transport 



The equations of ideal magnetohydrodynamics (MHD) can be written in con- 
servative form as 



| + V.(pv) = 


(1) 


dpV + V- (pvv BB) + VP* = 

at 


(2) 


<9B 

— + Vx(B x v)=0 
at 


(3) 


+ V • {{E + P*)v - B(B • v)) = 


(4) 
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where p is the mass density, pv the momentum density, B the magnetic field, 
and E the total energy density. The total pressure P* = P + (B • B)/2 where 
P is the gas pressure. This system of equations is closed with the addition 
of an equation of state which relates the pressure and density to the internal 



Throughout this paper we will assume an ideal gas equation of state for which 
P = (7 — l)e, where 7 is the ratio of specific heats. Note that we have chosen 
a system of units in which the magnetic permeability \x = 1. 

In addition to the evolutionary conservation laws, equations (1) through (4), 
the magnetic field must also obey the divergence free constraint, i.e. V • B = 0. 
Here, as in [1], this is accomplished using the method of constrained transport 
(CT). In this method one starts from the differential form of the induction 
equation (4) and constructs an integral relation by area averaging the compo- 
nents of B normal to the grid cell faces over the respective face and applying 
Stoke's theorem. Analogous to the finite volume method, the resulting integral 
relation forms the basis of the numerical evolutionary equation. Note that one 
immediate consequence of CT is that the V • B = constraint on the magnetic 
field is satisfied in an integral sense over the smallest discretization scale, the 
grid cell. The volume averaged magnetic field components, which for example 
are necessary to calculate the internal energy, are defined equal to the average 
of the interface averaged components. 

In this paper we will assume a regular, three dimensional, Cartesian grid. We 
will use the standard notation that grid cell (i,j,k) is centered at (xi,yj,Zk) 
and has a size (5x,5y,5z). Time levels will be denoted by a superscript and 
interface values will be denoted by half increments to the index, e.g. the volume 
averaged x-component of the magnetic field at time t n is defined to be 



3 Calculating the Interface States 

In this section we describe the calculation of the "interface states" in the PPM 
algorithm for ideal MHD in three dimensions. The PPM interface state algo- 
rithm is based upon the idea of dimensional splitting, and as a result is a 
one-dimensional algorithm including both spatial reconstruction and a char- 
acteristic evolution of the linearized system in primitive variables. For ideal 
MHD, however, it was shown in [1] that it is necessary to include multidi- 



energy, 



e = £-p(v-v)/2-(B-B)/2 . 



(5) 




(6) 
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mensional terms when calculating the interface states. The three dimensional 
interface state algorithm is thus a generalization of the two-dimensional al- 
gorithm which for consistency must reduce to the two- and one-dimensional 
algorithm in the appropriate limits. The interface states in the PPM algorithm 
are typically calculated by evolving the system of equations in primitive vari- 
ables. Consider the induction equation, which in component form is 

dB x d d 

~dT + dy ^ VyBx ~~ ByVx > + ~d~ z ( VzBx ~ BzVx > = ^ 
dB„ d , d . 

~dt + dx~ {VxBy ~ BxVy) + d~z {VzBy ~ BzVy) = (8) 

dB z d . d . 

^ + — (v x B z - B x v z ) + — (v y B z - B y v z ) = . (9) 



In these equations there are terms proportional to dB x /dx, dB y /dy and 
dB z /dz which we will refer to as "MHD source terms". (When the system of 
equations for MHD is written in primitive variables, these source terms only 
appear in the induction equation. As a result we will not discuss the remain- 
ing MHD equations in this section.) The question before us is: which terms 
in the induction equation need to be included in the calculation of the inter- 
face states? In what follows we specialize to the calculation of the ^-interface 
states; the y- and ^-interface state calculation follows by symmetry. 



3.1 2D MHD Interface State Algorithm 



Before constructing the three-dimensional interface algorithm, it is instruc- 
tive to recall the two-dimensional algorithm presented in [1]. For the two- 
dimensional (x, y)-case, the induction equation for B z was simplified using the 
V ■ B = condition eliminating the MHD source terms from the evolutionary 
equation for B z . The result is the following set of equations for calculating the 
rc-interface states in 2D 

9B X , x 

-af=° (10) 

-gf + g^ M. " B X V V ) = (11) 

dB, d , , dv z , , 

_£ + _ feB J_ Bi _£ = . (!2) 



Utilizing the V • B = condition to eliminate the source terms from the evolu- 
tionary equation for B z bestows a very important property on the calculation 
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of the interface states, namely the balance of the MHD source terms is exactly 
and explicitly included. The importance of this step can be easily understood 
by considering the advection of a magnetic field loop initially confined to the 
(x, y)-plane (i.e. B z = 0) with v = constant, v z ^ 0, and (3 = 2P/B 2 3> 1. 
If the MHD source terms had not been eliminated from equation (9) prior to 
dimensional splitting and the following equation 

OB 

^ + — (v x B z - B x v z ) = (13) 

was used instead of equation (12) the ^-interface state would include an erro- 
neous B z evolution owing to the term v z (dB x /dx). Experience shows that this 
error is not eliminated when updating the interface states due to transverse 
flux gradients in the CTU algorithm, leading to steady, secular growth of B z 
which effectively warps the field loop. The important point to note here is that 
the balance of the MHD source terms resulting from the V ■ B = condition 
must be accurately represented in the calculation of the interface states. 



3.2 3D MHD Interface State Algorithm 



The 3D interface algorithm we construct is designed to explicitly incorporate 
the potential balance between the MHD source terms and to reduce exactly 
to the 2D interface states algorithm in the limit that the problems is two- 
dimensional and grid aligned. The essential idea is to rewrite the induction 
equation as follows prior to applying the idea of directional splitting. 



9B X f d . „ „ . 
~W + \~d~y ^ x ~ yVx) ~ Vx xy 
d_ 

dz 




+ I (v z B x - B z v x ) - v x L xz [^-v\\ = (14) 



dB { d 

~df + )dx ^ VxBy ~ BxVy) ~ VyLyx 
d_ 

dz 




+ 1 iL ( v * B v - B » v v) - v v L v* ( ) \ = ( 15 ) 



9B Z [ d . „ „ . fdB v 



6 



(16) 



where we've added a limited amount of the transverse MHD source term to 
each component of the electric field gradient and grouped terms according to 
the fashion in which they will be split. The mathematical form of the limiter 
functions, e.g. L xy , is determined by imposing constraints on the directionally 
split and unsplit system. Clearly, to recover the induction equation we have 



et cetra. Directionally split, we obtain the following system for the x-coordinate 



direction 






dB x 
dt 


= 




dBy 

dt 


+ ^ (v x B y 


- B x Vy) 


dB z 
dt 


+ §- x (v x B z 


- B x v z ) 



dB z 



->yx 



dz 
dB 



(19) 



dy 



v - = . (20) 



To determine the form of the limiter functions, we minimize the magnitude of 
the sum of the MHD source terms. For equation (20) we find 

L Jf)"'(-f>f) 



where the minmod function is defined as 



minmod (x, y) 



sign(:r) min(|x|, \y\) if xy > 
otherwise. 



(22) 



Note that this limiter function satisfies the constraint identified in equation 
(17). The mathematical form of the remaining limiter functions in equations 
(14 - 16) is given by cyclic permutation of (x, y, z) in equation (21) and appli- 
cation of the constraint noted in equation (17). 

There is also a simple physical argument for why the limiter function takes 
the form described by equation (21). Considering equation (9), if (dB x /dx) 
and (dBy/dy) have opposite signs, but not necessarily the same magnitude 
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we wish to incorporate the balance of these two MHD source terms by adding 
and subtracting the term with the smaller magnitude so that the resulting 
(reduced) MHD source term is associated with only one of the flux gradients. 
If, on the other hand, these derivatives have the same sign, then there is 
no balance between the source terms and the induction equation should be 
unmodified. This is precisely the result of the minmod limited source term in 
equation (21). 

Finally, we note that using the properties of the minmod function and the 
V • B = condition, equations (18-20) can be simplified to 

dB x , . 

-er = ° (23) 

9B y d , „ s „ dv v . , / dB x dB v \ „ 

+ ^ (v x B y ) n, -f d U±, U (24) 



dB * 9 < n\ n dv * ■ a ( dB * dB A n (on 

~w dx {VxBz) - Bx ^ ~ Vzmmmod Up —Br r (25) 



for calculating the x-interface states. As a practical matter, these limited MHD 
source terms are evaluated in terms of the cell average of the magnetic field 
gradients, i.e. for equation (25) in cell (i,j,k) we use 

minmod ( Bx > i+1 / 2 ' j > k ~ ^-V 2 ^ ^,j,fc-i/2 - B x ^ j>k+1 / 2 \ 
\ Sx 5z I 



The equations for the y- and z-interface states follow from cyclic permutations 
of (x, y, z). In the limiting two dimensional case of either d/dy — or d/dz — 
this approach reduces to the interface state algorithm outlined in §3.1. More- 
over, in the limiting two dimensional case of d/dx = 0, the x-interface state 
will equal the cell center state, just what one expects from one and two di- 
mensional calculations. As a result, this algorithm for calculating the interface 
states preserves the B z = condition for the gedanken experiment described 
in §3.1 involving the advection of a magnetic field loop. 



4 Constrained Transport (CT) Algorithm 



The coupling of a Godunov, finite volume algorithm with the method of CT 
requires an algorithm for constructing the grid cell edge averaged electric fields 
(emfs) from the Godunov fluxes. This algorithm is typically referred to as a 
CT algorithm. The process of applying a CT algorithm to calculate the CT 
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emfs from the Godunov emf's can be described as a predictor / corrector 
process where the Godunov emf's are the predictor values and the resulting 
CT emfs are the corrector values. In [1] a simple framework for constructing 
CT algorithms was presented and a few CT algorithms were constructed and 
tested with the 8 C CT algorithm selected as having the best properties. The 
£ c CT algorithm is constructed in such a way as contain an upwind bias 
(according to the contact mode) and to reduce to the correct Godunov emf 
for grid-aligned, plane-parallel flows. In this paper we will also use the S c CT 
algorithm which for the sake of completeness we briefly review here. 

Consider for the moment the calculation of the ^-component of the electric 
field at the grid cell edge (i + 1/2, j + 1/2, k). (The calculation of the x- and y- 
components of the CT electric fields follows an analogous procedure.) The CT 
algorithms described in [1] compute the CT electric field at this location from 
the four neighboring face center electric field components (Godunov fluxes) as 
well as estimates of the gradients of the electric field as follows 



£ Z ,i+l/2,j+l/2,k = ^ (£z,i+l/2,j,k + £z,i+l/2,j+l,k + £z,i,j+l/2,k + £z,i+l,j+l/2,k) 



+ Sy_( ( d£ 



Sx I ( d£ 



dx 



i+l/2J+l/4,fc 



j+l/4j+l/2,fc 



dy 

'd£z_ 
dx 



i + l/2,j + 3/4,k; 



i+3/4J+l/2,fc, 



(27) 



To complete this CT algorithm we need to specify a way to calculate the 
derivatives of S z on the grid cell face. The S c CT algorithm computes the 
electric field gradient at the grid cell face by selecting the "upwind" direction 
according to the contact mode, i.e. 



'd£z 
s dy 



i+l/2J+l/4,fc 



(d£ z /dy) i>j+1 / 4;k for u X)i+ i/ 2 j,fc > 

(d£ z /dy) i+ld+1 /4,k for v x>i+1 / 2 ,j,k < (28) 

j+i/4,k + (d£ z /dy) i+ljj+1 /4,k) otherwise 



with an analogous expression for the (d£ z /dx). The final detail involves the 
definition of the electric field derivatives in equation (28). These are com- 
puted using the face centered electric fields (Godunov fluxes) and a cell center 
"reference" value £ r z ^j : ki e -S- 

( =2 ( £*jJ+W ~ £z,ij,k \ / 2g x 

\ d v )i, j+ i/4,k \ s v J ' 
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In the 2D MHD CTU algorithm described in [1] and the 3D version described 
here, the cell center reference electric field £ r zi jk is computed using the cell 
center state at an appropriate time level. For the first interface flux calcu- 
lation, using the interface states described in §3.2, the cell center reference 
electric field is computed using the cell center state qfj k when integrating 
from time t n to t n+1 . The calculation of the cell center reference electric field 
in subsequent steps of the integration algorithm use time advanced states and 
will be described later in connection with the integration algorithm. 



5 Corner Transport Upwind Algorithm 

In this section we are interested in applying the CTU algorithm to the system 
of equations for ideal MHD. The CTU algorithm was originally described 
by Colella [7] as an unsplit, 2D finite volume algorithm for solving hyperbolic 
systems of conservation laws. The 3D generalization of the CTU algorithm was 
subsequently presented by Saltzman [18]. The CTU algorithm is generally set 
within a predictor-corrector formalism and utilizes PPM [8] when applied to 
Euler's equations - the archetypical system. 

Prior to delving into the details of the numerical algorithms, it is worth while 
pointing out that the system of equations for ideal MHD differs from Euler's 
equations in non-trivial ways which are very important when applying the 
CTU algorithm to MHD. First, a straight forward application of the direc- 
tional splitting technique to the MHD equations in primitive and conservative 
variables results in an incompatible set of equations. This results in the need 
to incorporate source terms in the transverse flux gradient corrections in the 
2D and 3D CTU integration algorithm [1] since the PPM interface state algo- 
rithm uses the primitive variable form of the equations. Second, the treatment 
of the multidimensional MHD source terms in the MHD PPM interface states 
algorithm, described in §3, results in the need to incorporate source terms 
in the transverse flux gradient updates to the transverse components of the 
magnetic field at the interfaces. These details follow from the balance between 
multidimensional flux gradients imposed by the V • B = condition. 

In this section we present two variants of the CTU integration algorithm. In 
§5.1 we present a brief, functional description of the 3D CTU algorithm as 
described by Saltzman [18] which we refer to here as the 12-solve algorithm 
since it requires 12 solutions to the Riemann problem per zone per time step. 
We will discuss the challenges associated with adapting this algorithm to the 
equations of ideal MHD. However, for a variety of reasons, the principle one 
being the complexity of the algorithm, we will not present the algorithmic 
elements for the 12-solve MHD CTU algorithm in detail. In §5.2 we present a 
simple variant on the CTU algorithm which requires only 6 solutions to the 
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Riemann problem per zone per time step and describe this 6-solve algorithm 
in detail. We summarize with a discussion of the strengths and weaknesses of 
this algorithm relative to the 12-solve CTU algorithm as a prelude to §6 where 
we present a variety of results comparing the 6-solve and 12-solve MHD CTU 
algorithms. 



5.1 12-solve CTU 



In this subsection we present a functional description of the 12-solve CTU 
algorithm as constructed for Euler's equations. This serves the goal of making 
the discussion more self-contained as well as allowing us to directly point out 
where particular elements of the integration algorithm pose challenges when 
applied to ideal MHD. For a more detailed description of the algorithm, or 
the theoretical underpinnings, see [7,15,18]. 

We begin by choosing a numerical flux function T{q Ll q R ) which is assumed 
to return a suitably accurate solution for the flux obtained by solving the 
Riemann problem associated with and q R , the left and right states. The 
12-solve CTU algorithm can then be described as follows. 

Step 1, calculate the left and right PPM interface states q*L X ,i+i/2,j,ki Q.*Rx,i+i/2,j,ki 

qly,i, j+ i/2,k, Q*ny,ij+i/2,k> 9L,ij,fc+i/2 and QjRz,ij,k+i/2 and the associated interface 
fluxes 



F x,i+l/2,j,k = -Fx{<lLx,i+l/2,j,k-> 9 Rx, 2, j,k) (30) 
F y,i,j+l/2,k = Fy(QLy,i,j+l/2,ki Q*Ry,i,j+l/2,k) ( 31 ) 
F z,i,j,k+l/2 = Fz{<lLz,i,j,k+l/2-> lRz,i,j,k+l/2)- ( 32 ) 

Step 2, for each interface state calculate two interface states evolved by St/3 
of a single transverse flux gradient, i.e. 

^Lx,i+l/2J,k = 1*Lx,i+l/2,j,k + 3^ ( F y,i,j-l/2,k ~ F y,i,j+l/2,k) ( 33 ) 

9j&,i+l/2J,fc = QRx,i+l/2J,k + 3^ ( F y,i+l,j-l/2,k ~ F y,i+l,j+l/2,k) ( 34 ) 

lLx,i+l/2,j,k = Qlx,i+l/2,j,k + ^ ( F z,i,j,k-l/2 - F z,i,j,k+l/2) ( 35 ) 

QRx,i+l/2,j,k = Q.Rx,i+l/2,j,k + 3^ ( F z,i+l,j,k-l/2 ~ F z ,j ,k+l / 2) ( 36 ) 

with y- and z-interface states being defined in an equivalent manner by cyclic 
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permutation of (x,y,z) and (i,j,k). For each of the St/3 updated interface 
states, calculate the associated flux, giving the two x-interface fluxes 



p*\y _ <r (*\y *\y \ (n 7 \ 

r x,i+l/2,j,k ~ • r x\HLx,i+l/2,j,ki L lRx,i+l/2,j,k) \° ' ) 

F x ]i + i/2,j,k = 3~x( ( lLx,i+l/2,j,k-> lRx,i+l/2,j,k) (^) 

and similar expressions for the y- and z-interface fluxes. 

Step 3, at each interface evolve the PPM interface states by St/2 of the trans- 
verse flux gradients, i.e. 



n+1/2 * I 01 ( p*\z _ p*|z \ /oq\ 

L lLx,i+l/2,j,k t dLx,i+l/2,j,k 2Sy V f>*>J'-V2,fc r y,i,j+l/2,k) l oy ) 

+ oxT {^z!li,k-l/2 ~ ^,^7^+1/2) (40) 



"7^ \ i I z,y,*;-l/2 ~~ - t ' 2 ,i,i,fc+l/2 

^ n+1 /2 * , ( F *\z _ p *\z \ /.-.x 

< iRx,i+l/2,j,k — 'iRx,i+l/2,j,k^ 2Sy V r y,i+l,j-l/2,k r y,i+l,j+l/2,k) \* L ) 

+ 2S~Z (^+ 1 .J. fc -V2 ~~ ^M+lJ.fc+l^) ( 42 ) 

with y- and z-interface states being defined in an equivalent manner by cyclic 

permutation of (x,y,z) and (i,j,k). For each of the St/2 updated interface 
states, calculate the associated flux, giving the x-interface flux 

pn+l/2 _ T / n+1/2 n+1/2 \ {A o\ 

r x,i+l/2,j,k — "> x\HLx,i+l/2,j,k^ L lRx,i+l/2,j,k) J 



and similar expressions for the y- and z-interface fluxes. 

Step 4, update the the conserved variables from time n to n + 1 via the fully 
corner coupled numerical fluxes 



n n+i _ n n ,2L( F n+1/2 - F n+1/2 ) (U) 

%j,k — %j,k ^ fa \ r x,i-l/2,j,k r x,i+l/2,j,k) 

St_ I n+1/2 _ pn+1/2 \ St_ / „+i/ 2 ^ p n+l/2 \ /^x 

Sy V J/.^'-Va.fe r y,i,j+i/2,k) ^ ^2,^^-1/2 r z,i,j,k+i/V - 

This completes the description of the 12-solve CTU algorithm for a typical sys- 
tem of conservation laws, such as Euler's equations. Unfortunately, as written 
above, the 12-solve CTU algorithm does not result in a useful method for ideal 
MHD. This can be understood on rather general grounds by noting that the 
intermediate steps in the algorithm use partial updates based on a dimensional 
splitting of the equations in conservation form. This in turn ignores the poten- 
tial balance between flux gradients in different directions (in particular MHD 
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source terms associated with those flux gradients) which is always present for 
MHD owing to the V • B = constraint. We make this point more concrete 
by considering two points in detail. 

First, note that the parallel flux gradient terms (x-flux gradient at x-interfaces, 
etc.) are included in the PPM interface states algorithm using the dimension- 
ally split, primitive form of the equations for MHD. Meanwhile, the transverse 
flux gradient terms are included using the conservative form of the equations. 
Since the dimensionally split primitive and conservative form of the equations 
for MHD are not commensurate, this amounts to neglecting certain MHD 
source terms resulting in a formally first order accurate integration algorithm. 
In addition, such an algorithm would also show secular evolution of a mag- 
netic field component perpendicular to the magnetic field loop in the gedanken 
experiment discussed in §3. Note that the same is true of the two-dimensional 
CTU algorithm [1], which required the addition of source terms to the trans- 
verse flux gradient correction step. For the 3D 12-solve CTU algorithm, with 
two predictor steps, this source term correction procedure is increasingly com- 
plicated. 

To see why this is so, consider step 2 in the description of the 12-solve CTU 
algorithm just presented. In particular, note that in this step one generates 
two interface states at each interface, by evolving the PPM interface state by 
St/ 3 of one transverse flux gradient. In other words, for each interface normal 
component of the magnetic field (used to define the divergence of the mag- 
netic field), one generates two normal components, each of which is evolved 
by one half of a CT or Stokes EMF field loop. (As an aside, note that at 
this stage V • B = is satisfied in the sense of the average of these normal 
magnetic field components.) It is useful at this stage to consider again the 
gedanken experiment discussed in §3 and consider the generation of the two 
z-interface states in step 2 of the integration algorithm. As a result of the 
MHD source terms in the x- and y-flux gradients, namely v z (dB x /dx) and 
v z (dBy/dy), these two z-interface states will show a non-zero, in fact equal 
and opposite, evolution of B z . This is a manifestation of a failure to satisfy the 
balance condition discussed in §3. Thus what we find is that it is the balanced, 
dimensionally split system presented in §3, equations (14)-(16), which should 
replace the straight forward dimensionally split conservative form of the in- 
duction equation in steps 2 and 3 of the 12-solve CTU algorithm. In practice, 
this means that the interface normal components of the magnetic field must 
incorporate a source term so as to maintain this balance, and the two, say 
z-interface states, which are generated must use the same source term, with 
opposite sign, so as to maintain the magnetic divergence condition in an aver- 
age sense. Finally, note that as a result of dimensionally splitting the system, 
the momentum and energy update in step 2 and 3 also require the addition of 
source terms to balance terms like B x (dB x /dx) and B • v (dB x /dx) etc. 
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The net result is that a well balanced, 12-solve CTU algorithm for ideal MHD 
can be constructed by using partial updates based on a dimensional splitting 
of the MHD equations using a carefully chosen, non-conservative form for the 
intermediate steps. This form is found by paying particular attention to the 
implicit balance between the flux gradients, as was done for the induction 
equation in §3 for the PPM interface state algorithm. The advantage of this 
approach is a computational algorithm which is optimally stable for CFL 
numbers < 1. The disadvantage is that the algorithm is complicated. We 
have implemented the 12-solve MHD CTU algorithm as described above and 
present results of tests of the method in §6. However, the complexity of the 
method motivates us to find a simpler alternative, which we describe below. 



5.2 6-solve CTU variant for MHD 



In this subsection we present a simple variant on the 12-solve CTU algorithm 
which we will henceforth refer to as the 6-solve algorithm. For Euler's equa- 
tions, the 6-solve algorithm can be described concisely as the 12-solve CTU 
algorithm of §5.1 omitting step 2 and replacing F^ +1 , 2 - k and F*^ +1 , 2 ■ k with 
Fx,i+i/2,j,k ( an d similarly for the y- and ^-fluxes) in step 3. Alternatively, one 
may also describe it as a formal extension of the 2D CTU algorithm in which 
the parallel and transverse flux gradients are included in the interface states 
in a two-step process. In what follows we present a functional description of 
this 6-solve CTU algorithm for MHD including a detailed description of the 
treatment of the MHD source terms and constrained transport electric fields. 

Step 1, calculate the left and right PPM interface states Q*L X ,i+i/2,j,ki QRx,i+i/2,j,ki 
Qly,ij+i/2,k> q*Ry,ij+i/2,k> Q*Lz,i,j,k+i/2 and q Rz ^ k+1/2 including the MHD source 
terms described in §3.2 and the associated interface fluxes 



F x ,i+l/2,j,k = v(.Q.Lx,i+l/2,j,ki QRx,i+l/2,j,k) 
Fy,i,j+l/2,k = J~ y{ ( lLy,i,j+l/2,ki Q.Ry,i,j+l/2,k) 
Fz,i,j,k+l/2 = 3~'z(qLz,i,j,k+l/2i QRz,i,j,k+l/2) ■ 



(46) 
(47) 
(48) 



Step 2, apply the CT algorithm of §4 to calculate the CT electric fields 

^,ij+i/2,fc+i/2. £y,i+i/2,j,k+i/2 and £**i+i/2j+i/2,k usin S the numerical fluxes from 
step 1 and a cell center reference electric field calculated using the initial data 
at time level n, i.e. q^j k- 

Step 3, at each interface evolve the PPM interface states by 8t/2 of the trans- 
verse flux gradients. The hydrodynamic variables (mass, momentum and en- 
ergy density) are advanced using 
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ra+1/2 _ * | ^ f TP* ^ p* "\ 

( lLx,i+l/2,j,k — t iLx,i+l/2,j,k ' \ r y,i,j~l/2,k r y,i,j+l/2,k J 

+ 25z (^j*- 1 / 2 ~~ F li,j,k+i/2) + Y s x,i,j,k (49) 

n+l/2 _ * >lLfF* _ F* "\ 

lRx,i+l/2,j,k - L lRx,i+l/2,j,k ' \ r y,i+l,j-l/2,k r y,i+l,j+l/2,k) 

+ 26z { F *> i + 1 >i> k - 1 / 2 ~ F li+lJ,k+i/^) + ~2~ S x,i+l,j,k (50) 

where the x-interface MHD source term for the momentum density 

(S x ,iJ,k)pv = Bij,fe f J (51) 



dx , , 

/ i,j,k 



and the energy density 



(S Xtidtk ) E = (B y v y ) i:jjk mmmod f ^ j + 
, v . J dB v dB x \ 

{B z v z )ij :k mmmod I — , I . (52) 

The magnetic field components are evolved using the CT electric fields in place 
of the predictor fluxes. The interface normal component of the magnetic field 
is evolved using the integral form of the Stokes loop, 



R n+l/2 _ Rn (c* _ c* 

n x,i+l/2,j,k~ n x,i+l/2,j,k 2§y \°z,i+l/2,j+l/2,k °z,i+l/2,j-l/2,k 

+ 25z (^+V2J,fc+l/2 ~~ £y,i+l/2,j,k-l/2) ■ ( 53 ) 

The y-component of the magnetic field is evolved using 



(By) Lx/+l/2,j,k ~ (By)Lx,i+l/2,j,k ^ {^x,i,j+l/2,k+l/2 £x,ij+l/2,k-l/2) 

-—(f* -F* ) 

46 Z \ xjj-VW+VZ C x,i,j-l/2,k-l/2j 

+j(Sx, i , j , k ) By (54) 

(By) Rx/+l/2j,k = (Bv)lix,i+l/2,j,k ~ ~^£~ z {fx,i+\,j+l/2,k+l/2 ~ ^x,i+\,j+l/2,k-l/2) 
45Z V C ^+ 1 .J- 1 /2,fe+l/2 C x,i+l,j-l/2,k-l/2j 

+-?r(Sx,i+ij,k)B y (55) 
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with 




(S x ,i,j,k) By = (v y )i^ k mmmod [ - — , — ] . (56) 

i,j,k 



The z-component of the magnetic field is evolved using 



( B z)llJ+l/2,j,k = ( B z)*Lx,i+l/2,j,k + (^x,i,j+l/2,fe+l/2 ~ ^x,i,j-l/2,fc+l/2) 

+ 45y (^.*.J'+V2,fc-l/2 ~~ £x,M-l/2,fc-l/2) 

+ y(5'x,i,j,fe)B z (57) 
( B z)l^J + l/2,j,k = ( B z)*Rx,i+l/2,j,k + (^x,i+l,j+l/2,fe+l/2 _ £x,i+l,j-l/2,k+l/2) 

+ A6y (^+ 1 'J+ 1 /2,fe-l/2 ~ ^x,i+l,j-l/2,fc-l/2) 

+7r(^x,i+i,j,fc)s z (58) 



with 



{S x ,i,j, k )B z = (^)ij,feminmod > "^r^ ■ (59) 



Note that the origin of these MHD source terms for the transverse components 
of the magnetic field can be clearly seen as resulting from the directional 
splitting of the induction equation described in §3.2. The momentum and 
energy density MHD source terms originate from the use of the primitive 
variable form of the MHD equations to calculate the PPM interface states. 
The y- and z-interface states are advanced in an equivalent manner by cyclic 
permutation of (x,y,z) and (i,j,k) in the above expressions. 

Step 4, for each of the St/2 updated interface states, calculate the associated 
flux, giving the x-interface flux 

pn+l/2 _ T / Ti+1/2 n+1/2 x / fif .x 

r x,i+l/2,j,k - r x\HLx,i+l/2,j,ki l iRx,i+l/2,j,k) v UL V 



and similar expressions for the y- and z-interface fluxes. 

Step 5, apply the CT algorithm of §4 to calculate the CT electric fields 

CS+ 2 i/2,ifc+i/2> 1 i/2J,fe+i/2 and OvWiAfc usin § the numerical fluxes from 
step 4 and a cell center reference electric field calculated using the cell aver- 
age state at time level n + 1/2 which is calculated as follows. The cell center 
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magnetic field components are defined as equaling the arithmetic average of 
the interface magnetic field components, B^k = (i^A,- fc + B ^W*)/ 2 
and similarly for the y- and ^-components. The mass and momentum density 
are computed using a conservative update with the predictor fluxes from step 
1, 



"+1/2 
%j,k 



+ 



St 



+ 



6t 



( T?* P* 

\ r x,i~l/2,j,k r x,i+l/2 



25y 



R, 



5t 



(61) 



y,i,j-l/2,k 



F y,i,j+l/2,k) + { F z,i,j,k-l/2 F z,i,j,k+^) 



Step 6, update the solution from time level n to n + 1. The hydrodynamic vari- 
ables (mass, momentum and energy density) are advanced using the standard 
the flux integral relation, 



n+l_n j_^Lf^ n + 1 / 2 _ pn+1/2 \ / fi9 x 

%j,k — %j,k -r y r x ,i-i/2,j,k r x,i+i/2,j,k) y ^) 

5t_ / n+1/2 pn+1/2 \ &t_ ( p n+l/2 n +l/2 \ 

<fy V y^j-V^ k r y,i,j+l/2,kJ ^2,^^-1/2 r z,i,j,k+l/2) 

and the interface averaged normal components of the magnetic field are ad- 
vanced using a Stokes loop integral, 



R n+1 

XJ x,j+l/2,i,fc 



^ t pn+1/2 



,71+1/2 



x,i+l/2J,fe r ^ 2 ,i+l/2j+l/2,fe 2,i+l/2j-l/2,fe) 
/„n-4-1/2 ~n+l/2 



02 V J/,»+l/2J,fe+l/2 °J/,i+l/2,j 



(63) 



R n+1 _ R n , (pn+1/2 

n y,i,j+l/2,k - n y,i,j+l/2,k ^ ^,j+l/2 j+l/2,fc 



C«+l/2 \ 
C 2 ,i-l/2,j+l/2,fcj 



_ ( pn +l/2 

fi z y C x,i,j+l/2,k+l/2 



pn+1/2 \ 
C x,i,j+l/2,k-l/2) > 



(64) 



and 



R n+1 _ R n ^ ( pn+1/2 _ p n+l/2 \ 

D z,i,i,k+l/2 — D z,ij,k+l/2 c y°y,i+l/2,j,k+l/2 °y,i-l/2,j,k+l/2 J 



+ 



Sx 

5t 



_ (pn+l/2 _ f n+l/2 \ 

§ y y°x,i,j+l/2,k+l/2 °x,i,j-l/2,k+l/2) ■ 



This completes the description of the 6-solve CTU algorithm. This relatively 
simple 3D integration algorithm is second order accurate and has the advan- 
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tage over the 12-solve CTU algorithm that no source terms need be included 
in the evolution of the interface normal components of the magnetic field. 
This algorithm is designed in such a way that for grid aligned flows it reduces 
exactly to the 2D CTU and ID PPM integration algorithms for problems in- 
volving the relevant symmetry. Additionally, consideration of the field loop 
advection gedanken experiment described in §3 shows that the 6-solve CTU 
algorithm is well balanced and preserves the B z = condition exactly. The 
downside of the present 6-solve algorithm is that we observe experimentally 
that the algorithm is stable for CFL < 1/2. When compared to the 12-solve 
algorithm, this is compensated by the fact that it requires half as many Rie- 
mann solutions per time step. Hence, to a large extent the 6-solve and 12-solve 
algorithms show similar computational cost: two time-steps with the 6-solve 
algorithm at a CFL number of 1/2 is nearly equivalent to one time-step with 
the 12-solve algorithm with a CFL number of one. 



6 Tests 

In this section we present results obtained with the 6-solve CTU + CT inte- 
gration algorithm just described. For the sake of comparison, and clarification 
of the dominant differences between the results obtained with the 6-solve and 
12-solve algorithms, some results using the 12-solve algorithm will also be in- 
cluded. We will find in this section, through a series of tests, that the dominant 
difference between the two are the stability domain. Otherwise, they are quite 
comparable in accuracy and computational cost. As a result, in practical ap- 
plications we prefer the 6-solve algorithm for MHD on account of its simplicity 
and smaller memory footprint. 



6. 1 Field Loop Advection 

In this section we discuss and present results for the advection of a magnetic 
field loop. In order to narrow the focus and clarify the discussion, in this section 
we will concern ourselves primarily with the the initial conditions in which the 
density p, velocity v, and pressure P are constants, and the magnetic field is 
weak in the sense that (3 = 2P/B 2 ^> 1. In this limit, the evolution equations 
for the magnetic field are well approximated by the advection of a set of passive 
scalar functions, say the components of the magnetic vector potential. In the 
construction of the 2D CTU-CT algorithm [1] as well as the 3D 6-solve and 
12-solve algorithms presented here we have found that recovering the correct 
solution in this limiting case can be surprisingly difficult for conservative, finite 
volume algorithms applied to the ideal MHD equations. 



18 



As a concrete example of a situation in which this problem can be challenging, 
consider the field loop advection test problem studied in [1] and discussed as 
a gedanken experiment in §3. Specifically, consider a field loop confined to the 
(x, y)-plane, i.e. B z = 0, and a constant advection velocity field with v z ^ 0. 
If care is not taken to respect the balance between the MHD source terms in 
calculating the interface states, updating them with transverse flux gradients, 
etc. one can find an erroneous and sometimes secular evolution of B z . It should 
be noted that these concerns are not limited to CTU or PPM which use a 
predictor step. For example schemes with are conservative but do not satisfy 
V ■ B = will also find erroneous evolution for B z in this problem. There is 
an interesting corollary which results from this observation. If a conservative 
numerical algorithm can solve this magnetic field loop advection problem and 
preserve the solution B z = for all time, it also satisfies the V ■ B = 
condition. 

We begin by selecting a computational domain —0.5 < x < 0.5, —0.5 < y < 
0.5, and — 1 < z < 1, resolved on a N x N x 2N grid and apply periodic 
boundary conditions. The hydrodynamical state is uniform with a density 
p = 1, pressure P = 1, and velocity components (v x , v y , v z ) = (1, 1, 2). 
The initialization of the magnetic field is most easily described in terms of a 
vector potential in the coordinate system (x\, X2, £3) which is related to the 
computational coordinate system (x, y, z) via the rotation 

xi = (2x + z)/VE 
x 2 = y 

x 3 = (-x + 2z)/V5 . (66) 
In particular, we choose A\ = A 2 = and 



A, 



B (R-r) for r < R 

1 ; " (67) 
for r > R 



where B = 10~ 3 , R = 0.3 and r = ^Jx\ + x\ in the domain — 0.5Ai < x\ < 
0.5Ai, — 0.5A 2 < x 2 < 0.5A 2 . To satisfy the periodic boundary conditions 
we choose Ai = 2/v^5 and \ 2 = 1 and define A 3 (xi + n\\,x 2 + m\ 2 ,x 3 ) = 
A 3 (xi, x 2 , X3) for all integers (n, m). 

As a quantitative measure of the dissipation in the algorithm we plot the 
time evolution of the volume averaged magnetic energy density < B 2 > nor- 
malized to the initial (analytic) value < B 2 >= {\fhnR 2 /2)B 2 in Figure 1. 
Interestingly, the magnetic energy density < B 2 > shows a temporal evolu- 
tion similar to what was observed with the 2D algorithm in [1]. Namely, it 
can be well fit as a power law of the form < B 2 >= C(l — (t/r) a ) where 
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t = (3.22 x 10 2 , 3.68 x 10 3 , 2.65 x 10 4 ) and a = (0.365, 0.328, 0.320) for 
N = (32, 64, 128) respectively. Moreover, these values are quite comparable 
to the time constant r = 1.061 x 10 4 and exponent a = 0.291 found in the 2D 
calculation. 

For the specific case of a cylindrical magnetic field loop with translation invari- 
ance in the z-direction (d/dz = 0) the 6-solve and 12-solve algorithms studied 
in this paper reduce exactly to the 2D algorithm presented in [1]. As such, 
the 3D algorithms presented here give the same solution as the 2D algorithm 
in [1] and preserve the solution B z = for all time (we have explicitly tested 
that this is true with our implementation of the method). With the axis of 
the cylindrical field loop is aligned along a non-special direction with respect 
to the grid, preserving this property is non-trivial. 

As a quantitative measure of the ability of the algorithm to preserve B3 = we 
plot the normalized error < \B 3 \ > /B in Figure 1 . This error is calculated by 
contracting the cell center magnetic field with a unit vector in the x 3 -direction 
and computing the volume average of its absolute value. From this plot it is 
clear that the convergence rate of < \B%\ > /Bq as measured in either the 
initial conditions, or the solution at time = 1 is approximately first order. 
This behavior is consistent with the observation that the 3-component of the 
magnetic energy is dominant on the axis and at the boundary of the magnetic 
cylinder where the current density is initially singular, as shown in figure 2. It 
is also worth noting that away from these regions, the solution preserves the 
3-component of the magnetic energy quite small. This would not be the case 
if care were not taken to balance the MHD source terms in the integration 
algorithm. 



Normalized <B > Evolution Normalized < IB 3 I > Evolution 
~~ ' 1 1 1 1 1 1 1 O.Olr i 1 i 1 i 1 i 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

time time 



Fig. 1. Time evolution of the normalized, volume average magnetic energy density 
< B 2 > and the component along the ^-direction, < B 2 >, for three different grid 
resolutions using the 6-solve integration algorithm. 
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Fig. 2. Thresholded image of the magnetic energy (left) and the 3-component of the 
magnetic energy, -Bf/2, at time = 1. 

6.2 Linear Wave Convergence 

In this subsection we present the results of a convergence study for both the 
6-solve and 12-solve CTU-CT MHD algorithms. The problem we study is the 
propagation of linear amplitude, planar waves in a direction which is oblique to 
the grid. The physical conditions of the problem are most easily described in a 
coordinate system (xi, x 2 , x 3 ) which is chosen such that the wave propagates 
parallel to the Xi-axis. In this coordinate system, the initial conserved variable 
state vector is given by 



where q is the mean background state, e = 10 is the wave amplitude, and R p 
is the right eigenvector in conserved variables for wave mode p (calculated in 
the state q). In order to enable others to perform the same tests presented here 
and compare the results in a quantitative manner, we include the numerical 
values for the right eigenvectors in the appendix. 

The mean background state q is selected so that the wave speeds are well 
separated and there are no inherent symmetries in the magnetic field orien- 
tation (when initialized on the grid). The density p — 1 and gas pressure 
P = I/7 = 3/5. The velocity component v 1 = 1 for the entropy mode test 
and vi — for all other wave modes. The transverse velocity components 
v~2 — V3 = 0. The magnetic field components B 1 — 1, B 2 — 3/2, and B 3 = 0. 
With this choice, the slow mode speed c s = 1/2, the Alfven speed c a — 1, and 
the fast mode speed Cf = 2 in the Xi-direction. 




(68) 
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The computational domain extends from < x < 3.0, < y < 1.5, and 
< z < 1.5, is resolved on a 2N x N x N grid and uses periodic boundary 
conditions. Initializing this problem on the computational grid is accomplished 
by applying a coordinate transformation 



x = x\ cos a cos (3 — X2 sin (3 — x% sin a cos (3 
y = x\ cos a sin (3 + x 2 cos (3 — x 3 sin a sin /3 

z = X\ sin a + x 3 cos a (69) 

from the (x\, x 2 , x%) coordinate system to the (x, y, z) coordinate system 
of the grid with sin a = 2/3 and sin/5 = With this choice, there is one 

wave period along each grid direction and the wavelength A = 1. The interface 
components of the magnetic field are initialized via a magnetic vector potential 
so as to ensure V • B = 0. 



The error in the solution is calculated after propagating the wave for a distance 
equal to one wavelength. Hence, the initial state is evolved for a time t = 
A/c where c is the speed of the wave mode under consideration. For each 
component s of the conserved variable vector q we calculate the Ll-error with 
respect to the initial conditions 

S <ls = ^3 E - <4 M I ( 7 °) 



by summing over all grid cells (i,j,k). We use the cell center components of 
the magnetic field in computing this error. In figure 3 we plot the norm of this 
error vector 



INI = JE(^) 2 (7i) 



for the fast, Alfven, slow and entropy modes. Both algorithms demonstrate a 
second order convergence. With the exception of the slow mode, the 6-solve 
algorithm shows lower errors than the 12-solve algorithm. Note that the choice 
of maximum resolution in the convergence study for each algorithm and wave 
mode was selected on the basis of the "cost" of the computation. 



6.3 Circularly Polarized Alfven Wave 



In this section we present results for the propagation of a circularly polarized 
Alfven wave in a periodic domain using both the 6-solve and 12-solve algo- 
rithms. This problem is interesting from the perspective that the wave is an 
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Fig. 3. Linear wave convergence of fast, Alfven, slow and entropy modes using the 
CTU + CT 6-solve and 12-solve integration algorithms. The symbols denote the 
calculated Ll-error norm. 

exact nonlinear solution to the ideal MHD equations. Hence, this problem en- 
ables one to easily measure the nonlinear convergence to a multidimensional 
solution of the ideal MHD equations [20]. This problem is also interesting 
from the point of view that for a range of parameters, the circularly polarized 
Alfven wave is susceptible to a parametric instability [12,10]. Unfortunately, 
this situation has also hindered its applicability as a general and robust test 
for multidimensional MHD algorithms [16,14]. For the parameters used here, 
and suggested by Toth [20] we find no indication of instability. 

As with the linear wave propagation study presented in §6.2, the initial con- 
ditions are most easily described in a coordinate system (xi, x 2 , xs) which 
is chosen such that the wave propagates parallel to the xi-axis. In this coor- 
dinate system, the magnetic field components B 1 — 1, B 2 — 0.1 sin(27nri/A), 
and -B3 = 0.1 cos(27ra;i/A). The velocity components v\ = (0,1) for travel- 
ing or standing Alfven waves respectively, v 2 = 0.1 sin(27nci/A), and 1*3 = 
0.1 cos(27ra; 1 /A). The mass density p — 1 and the gas pressure P = 0.1, hence 
P = 2P/B 2 ~ 0.2. 

The computational domain used in this section is identical to that used in §6.2. 
In particular, we use the coordinate transformation given by equations 69 and 
a magnetic vector potential to initialize the magnetic fields so as to ensure 
V ■ B = 0. It is worth noting that this approach will necessarily result in mag- 
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netic pressure perturbations as a result of truncation error in initializing the 
magnetic field on the grid. The parallel component B\ is a constant, and hence 
rotation of this component will still result in a constant set of field compo- 
nents with no "pressure" variation. The perpendicular components (B 2 , B 3 ) 
however, will suffer some truncation error on initialization. Since B\jP = 0.1 
this truncation error in initialization will drive compressive waves. Note that 
with this set of initial conditions and v\ — the Alfven wave will travel a 
distance of one wavelength A in a time t — 1. 



As a quantitative measure of the solution accuracy, we present in figure 4 the 
norm of the LI error vector (as defined in equation 71) after propagating for a 
time t = 1 for both standing and traveling wave modes. From this figure we see 
that both traveling and standing circularly polarized Alfven waves converge 
with second order accuracy for both integration algorithms. The traveling 
wave mode shows a larger error amplitude relative to the standing mode, but 
it is worth noting that (while not shown here) the increase is fairly uniform 
over the components of the error vector. The 6-solve and 12-solve algorithms 
show quite comparable errors for both standing and traveling wave modes. 
When using a CFL number of 0.4, the 12-solve and 6-solve algorithms result 
in nearly identical errors. Increasing the CFL number to 0.8 with the 12-solve 
algorithm results in a slightly reduced traveling wave error, and increased 
standing wave error. These results indicate that the dominant difference in 
the LI error between the 6-solve and 12-solve algorithms results from the 
CFL dependence of the truncation error. 
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Fig. 4. LI error norm for the 6-solve and 12-solve integration algorithms for both 
standing (left) and traveling (right) circularly polarized Alfven waves. In particular 
note the dominant difference between the 12-solve and 6-solve errors is attributable 
to the CFL dependence. 

As a qualitative measure of the solution accuracy, we present in figure 5 scat- 
ter plots of B<i versus x\ for both standing and traveling wave modes after 
propagating for a time t — 5 using the 6-solve integration algorithm. These 
plots are constructed using the cell center components of the magnetic field, 
the cell center position and the coordinate transformation given by equations 
69. As a result of the fact that the wave is rotated with respect to the grid, 
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there are many grid cells with the same cell center xi-position. Hence, since 
these plots include every grid point in the grid, the lack of scatter in the plots 
demonstrates that the Alfven waves retain their planar symmetry throughout 
the calculation. Unfortunately, it is difficult to use the results presented here 
to make direct contact with solutions presented in the literature due to the 
scarcity of published three-dimensional test solutions. For analogous plots in 
a two-dimensional system see [20,16,1] 



Standing Alfven Wave Traveling Alfven Wave 




0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 



x i x i 

Fig. 5. Plot of B2 versus x± at t = 5 for the standing (left) and traveling (right) 
circularly polarized Alfven waves using the 6-solve integration algorithm. For com- 
parison, the initial conditions at t = for the N = 64 case is also included. 

As a final measure of the solution accuracy and convergence, we present re- 
sults for the dissipation of magnetic helicity in the case of a traveling circu- 
larly polarized Alfven wave. We note that this is not the cleanest possible 
test, since with periodic boundary conditions and a mean magnetic field, it 
does not appear to be generally possible to define a magnetic helicity which 
is conserved [5]. Nevertheless, we find that following [4] the magnetic helicity 
evolution associated with the fluctuating components of the magnetic field 
gives an interesting constraint on the problem considered here. In particular, 
let Bo =< B > (where angle brackets denote a spatial mean) and b = B — Bo 
denote the mean and fluctuating components of the magnetic field respectively. 
Also, define the magnetic vector potential associated with the fluctuating field 
as b = V x a. It is worth noting that as a result of periodic boundary con- 
ditions, B is time independent and the magnetic helicity associated with the 
fluctuating field H =< b ■ a > is gauge invariant. It follows that the time 
evolution of the magnetic helicity is given by 

4-<b-a>=-2<E-b> (72) 
dt 



where E is the electric field. Assuming ideal MHD, this equation can also be 
written as 

4 < b • a >= -2B - < v x b > . (73) 
at 
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From this expression it is clear for a circularly polarized Alfven wave the 
magnetic helicity should be conserved with < b • a >= Bj_/k. 

In figure 6 we present the time evolution of the normalized magnetic helicity 
H = (k/Bj_) < b a > for a traveling Alfven wave using the 6-solve integration 
algorithm for a variety of resolutions. The plots in this figure show two basic 
phenomena, dissipation and weak oscillations. The oscillations are an indica- 
tion that the circularly polarized Alfven wave is not resolved exactly. As such, 
certain features regarding the oscillations are worth mentioning. First, the os- 
cillation period r = 1/2 independent of the grid resolution and whether the 
Alfven wave is standing or traveling with respect to the grid. Second, the am- 
plitude of the oscillations in the helicity varies with resolution proportional to 
N~ 2 . Third, the oscillations are consistent in both amplitude and phase with 
the independently measured volume average quantity B - < v x b >. These 
details support the conclusion that the oscillations are a result of truncation 
error in resolving the circularly polarized Alfven wave. 

Normalized Helicity 
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Fig. 6. Plot of the normalized magnetic helicity H = (k/B 2 ^) < b-a > as a function 
of time for different resolutions. 

6.4 MHD Riemann Problem Inclined to the Grid 

In this section we present results for the solution of an MHD Riemann prob- 
lem in a three dimensional domain. The Riemann problem is a favorite test 
problem for computational algorithms since it can be chosen to study smooth 
flows, discontinuous flows, or a combination thereof. Moreover, the solution 
to this problem can, at least in principle, be calculated exactly allowing one 
to verify the algorithm in some parameter regime. For multidimensional algo- 
rithms this can also be an interesting test problem when the Riemann problem 
interface normal direction is chosen such that it has no special orientation with 
respect to the computational grid. In this configuration it provides a measure 
of the ability of the computational algorithm to faithfully reproduce the one- 
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dimensional solution on the large scale, despite the fact that on the scale of 
grid cells the flow contains multidimensional, interacting waves. In what fol- 
lows we describe the initial conditions, translation symmetry and boundary 
conditions for this problem and present the solution in a three-dimensional 
domain using the 6-solve CTU algorithm. 

We begin by choosing a coordinate system (x±,X2, £3) with the Riemann prob- 
lem interface located at X\ — and will use the terms left and right states 
to refer to the regions x\ < and x\ > respectively. To map the initial 
conditions to the computational domain, we apply the coordinate transforma- 
tion in equation 69 with the choice of rotation angles described below. This 
coordinate transformation can be inverted to read 



x 1 — x cos a cos (3 + y cos a sin (3 + z sin a 
£2 — —x sin (3 + y cos (3 

x 3 = —x sin a cos (3 — y sin a sin (3 + z cos a . (74) 

Using the fact that the initial conditions and solution to the Riemann problem 
are a function of the a;i-coordinate alone, the solution vector g(x + s) = g(x) 
for a translation vector s which satisfies xi(x + s) = xi(x). Making use of 
equation (74) we find that the continuous set of translation vectors s, for 
which the solution is invariant, satisfies the equation 

s x cos a cos [3 + s y cos a sin (3 + s z sin a = . (75) 



Now, for the problem at hand we are interested in the discrete set of trans- 
lation vectors for which (s x ,s y ,s z ) = (n x 5x,n y Sy,n z Sz) where {n x ,n y ,n z ) are 
integers and (5x, Sy, Sz) are the grid cell size in each direction. Making this 
substitution, and rearranging terms we find 

Sy Sz tan a . , 

n x + n v — tan/3 + n z - ^ = 0. (76) 

Sx Sxcosp 



We next choose the rotation angles (a, (3) such that 
Sy r x , Sztana r x 

/tan/5 = — and = — 77 

Sx r y Sx cos (3 r z 



where (r x ,r y ,r z ) are integers. With this choice, our equation for translation 
invariance becomes 

VjL + "hi + ^ = . (78) 

nn nr» nrt 

1 x 1 y 1 z 
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This is the key relation describing the discrete translation invariance of the 
initial conditions, and solution. 

There are a couple of interesting implications of this equation which are of 
practical importance for this Riemann problem. First, note that the transla- 
tion invariance described by equation (78) was constructed by considering a 
point translation symmetry and as such applies equally well to volume and 
interface averaged quantities. That is, there are no approximations involved 
in the statement that q id>k = qi+ nx ,j+n y ,k+n z for (n x ,n y ,n z ) which satisfy (78). 
Second, note that one coordinate direction, say the x-direction, can be iso- 
lated as the principle simulation direction and the transverse directions can 
be made as small as (r y ,r z ). Finally, note that the translation invariance re- 
lation (78) is the key relation for mapping computational grid cells to ghost 
cells for imposing boundary conditions. 

The specific Riemann problem we consider in this section is presented in [17] 
in test problem 2a. In the (x±, x 2 , x 3 ) coordinate system, the left state is 
initialized with p = 1.08, (v u v 2 , v 3 ) = (1.2, 0.01, 0.5), (B u B 2 , B 3 ) = 
(2/ v / 4tt, 3.6/Vtir, 2/ v / 47r) and P = 0.95. The right state is initialized with 
p = 1.0, (v u v 2 , v 3 ) = (0, 0, 0), (£?!, B 2 , B 3 ) = (2/v^F, 2/v^F) 
and P = 0.95. This problem is then mapped to the 3D domain with the 
rotation parameters (r x , r y , r z ) = (1, 2, 4). The computational grid has 
(Nx, Ny, Nz) = (768, 8, 8) grid cells covering the domain —0.75 < x < 0.75, 
< y < 1/64, < z < 1/64 and hence has a resolution of 5x = 5y = 5z = 
1/512. 

The solution to this Riemann problem at time = 0.2 is presented in figure 
7 using the 6-solve CTU algorithm. These plots include the cell-center data 
from every grid cell using the coordinate transformation in equation (69). The 
first thing to note in these plots is that since N y > r y and N z > r z there 
are multiple grid cells with the same cell-center xi-position. Therefore, the 
lack of scatter in these plots indicates that the algorithm retains the planar 
symmetry throughout the simulation. A comparison of the results presented 
here to the ID solution using the underlying PPM algorithm, with the same 
resolution, i.e. 5x = 1/512, indicates that the 3D solution has dissipation 
characteristics which are nearly identical to the ID algorithm. The dominant 
difference between the ID and 3D solutions is the presence of oscillations at 
the slow, Alfven and fast mode discontinuities. 

One question which has received a good deal of attention with this class of 
problem is the ability of the computational algorithm to maintain the par- 
allel component of the magnetic field, Bi, equal to a constant. We wish to 
point out here that oscillations are likely unavoidable unless the orientation 
of the Riemann problem is chosen to be aligned in a special direction with 
respect to the grid. As evidence of this fact, we note that the in the initial 
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Fig. 7. Solution to the Riemann problem in a direction oblique to the grid. 



conditions, the cell-center I^-component of the magnetic field shows an oscil- 
lation with an amplitude of approximately 8.26 x 1CT 3 despite the fact that 
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the interface averaged magnetic fields were initialized with an "exact" inte- 
gral average using a magnetic vector potential. This oscillation is therefore 
a result of the discretization relating the cell-center and interface averaged 
magnetic field components. In the initial conditions, as well as the solution at 
time = 0.2, the oscillations in B\ occur wherever the transverse components 
of the magnetic field rotate over a small scale such as the initial discontinuity, 
and the resultant fast, Alfven and slow mode discontinuities. Finally, we note 
that just as in the 2D paper [1], the oscillations in the parallel component 
of the magnetic field can be eliminated by restricting the solution to "macro- 
cells" . This operation effectively aligns the the Xi-direction with the macrocell 
[1, 1, 1] direction. 



6.5 MED Blast Wave 



Another problem which has been computed by a number of authors is the ex- 
plosion of a centrally over pressurized region into a low pressure, low (3 ambient 
medium. This is an interesting problem in the sense that it combines shocked 
flows, smooth flow regions, and strong magnetic fields. While the results of 
this test are not particularly quantitative in their measure of the accuracy, this 
test is a good measure of the robustness of the integration algorithm. Variants 
on this problem have been presented by a number of authors [21,3,13,1] and 
here we choose to use the parameters given by [13] for a three-dimensional 
domain. 

The computational domain extends from —0.5 < x < 0.5, —0.5 < y < 0.5 
and —0.5 < z < 0.5. The density p — 1, the velocity v = 0, and the magnetic 
field components B x = B z = 10/^/2 and B y = 0. Within a sphere of radius 
R = 0.125 about the origin the gas pressure P = 100 and (3 = 2P/B 2 = 2. 
Outside of this sphere, the gas pressure P = 1 and (3 = 2 x 10 -2 . These initial 
conditions are evolved until a time t = 0.02 using a 200 3 computational grid. 

In figure 8 we present images of the density, pressure, magnetic and kinetic 
energy density sliced along the y = plane at the end time. The general struc- 
ture of the solution is the same as one finds in the 2D calculation. Namely, 
the outermost surface in this expanding shell is a fast-shock which is only 
weakly compressive and energetically is dominated by the magnetic field. In- 
terior to this, one finds two dense shells of gas which propagate parallel to the 
magnetic field. These shells are bounded by a slow-mode shock and contact 
surface (separating the initially hot, interior gas from the surrounding cool 
ambient medium) on the outer and inner surfaces respectively. The maximum 
compression of the ambient gas by the slow-mode shock is approximately 3.3, 
the same as was found in the 2D calculation. The fact that the 2D and 3D cal- 
culations show quantitatively similar compression in the show-mode shock is 
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an indication that their motion is approximately one-dimensional, i.e. parallel 
to the magnetic field. 

The results of this section are interesting from the point of view that they 
demonstrate that the 6-solve integration algorithm is a robust algorithm, ca- 
pable of evolving shocked flows with (3 ~ 1CT 2 . Moreover, since the integration 
algorithm is unsplit, it preserves the symmetry of the initial conditions natu- 
rally. 




Fig. 8. Linearly scaled grey-scale images of the evolved state (time=0.02) for the 
MHD blast wave problem. The density (top left) ranges from 0.190 (white) - 2.98 
(black). The gas pressure (top right) ranges from 1.0 (white) - 42.4 (black). The 
magnetic energy density (bottom left) ranges from 25.2 (white) - 64.9 (black). The 
kinetic energy density (bottom right) ranges from 0.0 (white) - 33.1 (black). 
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7 Conclusion 



In this paper we have presented a three-dimensional MHD integration algo- 
rithm which combines the (6-solve) Corner Transport Upwind integration al- 
gorithm with the method of Constrained Transport for evolving the magnetic 
field. This integration algorithm is a natural extension, and generalization of 
the two-dimensional algorithm [1]. In addition we have outlined the essential 
elements to constructing a 12-solve CTU with CT integration algorithm for 
MHD and included results of this algorithm in §6. Both the 6-solve and 12- 
solve algorithms are found to be accurate and robust for approximately the 
same computational cost. As a result, we generally prefer the 6-solve algorithm 
as a result of its simplicity and smaller memory footprint. 

The three-dimensional MHD PPM interface states algorithm presented in this 
paper is a new and essential element of the integration algorithms. We have 
shown here that this is a natural extension of the 2D MHD PPM interface 
states algorithm presented in [1] and that it reduces identically to the 2D al- 
gorithm in the grid- aligned, plane-parallel limit. The 3D MHD PPM interface 
states algorithm was designed in such a way as to satisfy a multidimensional 
balance law involving what we have referred to here as MHD source terms. 
Failure to satisfy this balance law is found to result in erroneous and sec- 
ular evolution of the magnetic field under quite general conditions, e.g. the 
advection of a high (3 magnetic field loop. 

We have also presented a variety of test results for both the 6- and 12-solve 
MHD CTU CT integration algorithms. These test problems were selected so 
as to enable a comparison with previously published results, as well as to intro- 
duce new, quantitative measures of the the solution accuracy. One interesting 
result of these tests is the observation that the dominant difference in the LI 
error for the 6- and 12-solve algorithm convergence on smooth wave propaga- 
tion is attributable to the CFL number dependence. Throughout this section 
we have included the necessary information so as to enable other researchers 
involved in developing or applying MHD algorithms to make a quantitative, 
as well as qualitative, comparison with the results in this paper. 

Finally, it is worth noting that the integration algorithms presented here have 
been thoroughly tested on a great many test problems not included here. 
These include problems which are also of interest for their scientific merit. 
Examples include a study of the magneto-rotational instability [2] and the 
MHD Raleigh Taylor instability [19]. In a future paper we will detail our 
approach to combining the integration algorithms presented here with the 
methods of static and adaptive mesh refinement. 
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A Linear Wave Right Eigenvectors 



In order to enable others to perform the linear wave convergence test presented 
in section 6.2 and compare their results in a quantitative manner, we include 
the numerical values for the right eigenvectors here. In the wave-aligned coor- 
dinate system (xi, x 2 , x 3 ) the conserved variable vector and right eigenvectors 
(labeled according to their propagation velocity) are given by 
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